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With this work we present two new methods for the generation of thermostated, manifestly 
Hamiltonian dynamics and provide corresponding illustrations. The basis for this new class of 
thermostats are the peculiar thermodynamics as exhibited by logarithmic oscillators. These two 
schemes are best suited when applied to systems with a small number of degrees of freedom. 

PACS numbers: 



I. INTRODUCTION 

Back in 1984 Nose put forward a method for the gen- 
eration of equations of motion that sample the canonical 
ensemble pQ. The method is based on the Nose Hamil- 
tonian, reading: 

H = V Pl - + V(x) + + k B T\nX (1) 

i 

where a log-oscillator with Hamiltonian P 2 /2M + 
ksT In X, is non-linearly coupled to a "virtual" system 
(x,p). The thermostated dynamics of the "real" system 
are obtained after a time-rescaling and the application of 
a non-canonical transformation. The method was later 
further developed by Hoover [2], and is currently widely 
used and known as the Nose-Hoover thermostat. 

In this paper we unveil those special thermodynamic 
properties of log-oscillators which provide them with the 
power to act as thermostats and, based on them, show 
two more ways in which log-oscillators can be employed 
to generate thermostated dynamics. At variance with 
the method of Nose, these methods are genuinely Hamil- 
tonian, in the sense that the thermostated dynamics are 
obtained directly from Hamilton's equations of motion, 
with no need to perform a time rescaling nor the use of 
non-canonical transformations [2 Sj . Consequently these 
methods not only constitute a numerical means but, as 
well, can even be implemented in situ with real exper- 
iments aimed at thermostating a physical system. The 
first of the two methods has been reported recently with 
a letter, see in Ref. Its feasibility has been further 
discussed with a short account in Ref. HJ providing there 
the response which dispels a criticism raised by Hoover 
and co-workers [7]. 

It is important to stress that, just like the Nose-Hoover 
method, these methods only work provided the overall 
dynamics are ergodic, which might present a problem, 
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- especially when applied to small systems. In the case 
of Nose-Hoover thermostating one possible solution was 
offered by Martyna et al. [5J, who proposed the use of 
chains of Nose- Hoover thermostats. Our first method, 
at least in the implementation we have explored [that 
is considering a system of particles which interact with 
each other and with a log-oscillator via short range hard 
core repulsion, see [2l] below] seemingly is immune in ref- 
erence to this ergodicity issue |H1 HO] • Regarding our 
second method, see [29] below, the absence of ergodicity 
may present an issue; this second method, however, is 
sufficiently flexible as to overcome this challenge. 

II. HELMHOLTZ THEOREM 

The fact that logarithmic oscillators have a thermo- 
stating power is a consequence of their peculiar thermo- 
dynamic properties. In this section we shall clarify in 
what sense it is meaningful to talk about the thermo- 
dynamics of mechanical systems that have only one or 
few degrees of freedom, - as it is the case of logarith- 
mic oscillators, and demonstrate how to calculate their 
thermodynamic properties. 

Our starting point is the salient equation of thermo- 
dynamics: 

5Q/T = exact differential = dS (2) 

also known as the heat theorem [11] , As early as 1884, 
Helmholtz proved that this mathematical structure of 
thermodynamics is inherent to the classical Hamiltonian 
dynamics of systems having only one single trajectory 
for each energy, which he called monocyclic systems [12] . 
Arguably, this seldom appreciated and rarely known fact 
was one of the cornerstones on which ergodic theory 
(which generalizes Helmholtz monociclicity) and statis- 
tical mechanics were later built up by Boltzmann and 
others [nUMH]. 

The Helmholtz theorem goes as follows: Consider a 
classical particle in a confining potential <p(X; A), where A 
is an external parameter. To each couple (E, A) of values 



2 



of the energy and the external parameter is associated 
one closed trajectory in the system phase space. For 
each trajectory one can calculate the average quantities: 



*bT(S,A):=(£) 



F(E,X) 



<9A 



(3) 
(4) 



E,X 



where P, M are the particle momentum and mass respec- 
tively, and (•) e,x denotes time average over the trajectory 
specified by (E, A). Noticing that F(E, A) is the average 
force that the particle exerts against the external agent, 
keeping the parameter A at a fixed value, one realizes 
that 



SQ = dE + F(E, X)dX 



(5) 



represents the heat differential. The Hclmholtz theorem 
states that \/T(E 1 X) is an integrating factor for SQ, 



dE + F(E,X)dX 



T(E,X) 



exact differential = dS 



and that 
where 



S(E,X) = k B ]n<f>(E,X) 



(6) 



(7) 



2 / y/2M(E - <p(X; X))dX/h 

X_(E,X) 



= J dXdP6[E - H(X, P)]. 



(8) 



Here, X±(E, A) are the turning points of the trajectory, 
h is a constant with the units of an action, and 9{x) 
denotes Heaviside step function. Accordingly it is mean- 
ingful to call T(E, A) the temperature of the particle and 
S(E, A) its entropy. S(E, A) in [7] is also known as the 
Hertz entropy |15j . 

Once the function S(E, A) is known, one can then 
quickly calculate T(E,X) and F(E,X) in accordance to 
|6j as: 



/as 

[dE 



as /as 

~dX [dE 



(9) 
(10) 



and so obtain the thermodynamics of the system: equa- 
tion of state, specific heat, etc.. 

Following this scheme, in the next section, we will pro- 
ceed to derive the thermodynamics of log-oscillators and 
highlight the peculiar properties that provide them with 
thermostating power. 
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FIG. 1: (Color online) Black solid lines: Phase space trajec- 
tories of a log-oscillator at energies E = 1/2, 1, 3/2, . . . 9/2, 
inner curves have lower energies. Red dashed line: the mo- 
mentum distribution function, |17| Here M = 1, ksT = 1. 



III. THE PECULIAR THERMODYNAMICS OF 
A LOG-OSCILLATOR 

A. The heat capacity is infinite 

Let us consider a log-oscillator with Hamiltonian: 

H log (X,P) = ^+k B T\n 1 ^, (11) 

where M is the mass and b some positive constant with 
the dimension of length. Figure [T] depicts some trajec- 
tories in phase space of different energies. Solving the 
equation H\ og (X, P) = E for X, one sees that the trajec- 
tories are given by the equations: 



X = ± be E/kBT e -P>/ 



1Mk B T 



(12) 



That is the trajectories possess a Gaussian shape. Note 
that, accordingly, the maximal excursion grows exponen- 
tially with E/k B T: X max = be E/kBT . A straightforward 
calculation gives: 

$io g (£0 = J dXdP0[E - H log (X,P)} = 2b^2irMk B Te E/kBrj 

(13) 

Here, and in what follows we have set for convenience 
h = l. Accordingly, the entropy, [7J reads: 



S(E) 



E 
T 



k B ln[2b\/2TrMk B T}. 



Using the Hclmholtz theorem, we get: 

(P 2 /M) E = {dS/dE)- 1 = k B T. 



(14) 



(15) 



This expresses the major feature of the thermodynamics 
of a log-oscillator: all its trajectories inherit one and the 
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same absolute temperature, which is given by T, where 
T is the strength of the logarithmic potential. This fact 
is very peculiar: Consider for example the ID harmonic 
oscillator, in this case k B T(E) = E, namely the higher 
the energy, the higher the temperature. Similarly this is 
the case for a particle in a ID box, where kBT(E) = E/2. 

It therefore follows that the log-oscillator possesses a 
spectacular property: it has an infinite heat capacity; 



C(E) = (dT/dE)- 1 = oo 



(16) 



thus, it mimics a bath composed of an infinite collection 
of harmonic oscillators |16) , or one with an infinite num- 
ber of particles in a box. 



B. Log-oscillators sample the Maxwell distribution 

Yet another peculiar feature of the log-oscillator is that 
the probability density /(P) to find it with momentum 
P, is given by the Maxwell distribution at temperature 
T: 



f(P) = (2irMh B T) 



-1/2 -P 2 /2Mk B T 



(17) 



This holds independent of its energy E. To see this, con- 
sider the trajectory of the log-oscillator of some energy 
E. The probability to find the system at X, P, is given 
by the microcanonical distribution: 

P (X, P) = S[E - H log (X, P)]/n log (E) (18) 

where S(x) denotes Dirac's delta function, and 

n log (E) = J dXdP5[E-H log (X,P)} 

= dq> f^ E) = 2b^2^M/k B Te E ' kBT . (19) 

Therefore, the probability to find the log-oscillator at mo- 
mentum P, is obtained by the marginal distribution: 

/(P) = J dXp{X,P) = J dXS[E - H lo& (X,P)]/n los (E) 

=n^E)ml dx6 ^ E - p2 l 2M - kBT ' 



'In 



\x\ 



d 



&oxp[(E-P 2 /2M)k B T] 



fl log {E) 8E J 

2 d 



n log (E) dE 



be (E-P 2 /2M)/k B T = 



dX 

e -P 2 /2Mk B T 



y/2nMk B T 



(20) 



where we have used S(y) — d0{y)/dy. 

From 17 it is immediate to obtain that T(E) = 
(P 2 ) E /Mk B = T, in accordance with (15] 

The red dashed curve in Fig. [I] illustrates 17 When 
projecting the microcanonical distribution of the log- 
oscillator onto the P axis, the Maxwell distribution is 
obtained, regardless of the energy. 



IV. METHOD I 

The central feature of a thermal bath is that its heat 
capacity is infinite, hence, in this sense a single log- 
oscillator does indeed act like a thermal bath. Based 
on this observation it is reasonable to expect that when 
a system interacts weakly with a log-oscillator, the latter 
should induce thermostated dynamics at temperature T 
in the system. 

That this is indeed the case can be seen formally in the 
following manner [5]. Consider the total Hamiltonian: 

H (x, p, X, P) = H S ( X , p) + H log (X, P) + ft(x, X) (21) 

where 



ff s (x,p)=p 2 /2m + [/(x) 



(22) 



is the system Hamiltonian, and h(x, X) is a weak inter- 
action term that couples the system to the log-oscillator. 
Under the assumption that the total dynamics are er- 
godic, the probability density function p(x, p) for finding 
the system at (x, p) reads [17] : 



P(x,p) 



Hjogjgtot - ff s (x,p)] 

n(E tot ) 



(23) 



where E to t is the total energy of the compound system 
and 

fi(£tot) = J dXdPdxdp5{E tot -H(x,p,X,P)} (24) 

is the density of states of the compound system. Note 
that the shape of the distribution p(x, p) is given by the 
numerator, whereas the denominator only represents a 
normalization factor. Thus, from the fact that the den- 
sity of states of a log-oscillator is exponential in E/k B T, 



see 19 it immediately follows that: 



K X >P) 



3 -J? s (x,p)/fe B T 

W) 



(25) 



where Z{T) = /dxdpe- ffs(x ^ )/feBT . Thus, the constant 
temperature equations of motion read: 

x = p/m, 

p = -9 x (7(x) - cU(x, X) 
X = P/M 
[ P = -k B T/X -d x h( X ,X) 

where <9 X denotes the gradient operator in the x space and 
dx is a short notation for d/dX. Note that for h = 0, i.e., 
in absence of interaction, the system undergoes constant 
energy dynamics. 



Illustration 

Ref. [5] illustrates the numerical implementation of this 
method for small systems composed of few particles con- 
tained in a box and interacting through a repulsive hard 
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core potential 




(?r-(?)* 



\q\ > 2 1 / 6 (T 
■ E , |g| < 2 1 / 6 <7 



(26) 

The main limitation of this method comes from the fact 
that, in practical realizations, the logarithmic potential 
needs to be truncated at low values of A, for example by 
substituting it with: 



MX) 



k B T X- 
— - — In — 



b 2 



(27) 



This truncation results in a deviation of the single par- 
ticle velocity distribution from the target Maxwell dis- 
tribution. This deviation becomes more and more pro- 
nounced as the number of particles in the system in- 
creases, see Fig. 3 of Ref. El and can be compensated 
by rising the system energy as E ~ fkgT/2, where / 
is the number of degrees of freedom of the system. This 
energy rising, however, is accompanied by an exponential 
increase of the corresponding length and time scales in- 
volved in the dynamics which go as e E / klsT ~ e^ 2 , thus 
limiting the applicability of the method to systems with 
a small number of degrees of freedom. 

A prominent novel aspect of this method when com- 
pared to the other existing methods discussed in the lit- 
erature is that it can be implemented not only with com- 
puter simulations but also in analogue simulations, pro- 
vided one is able to implement the Hamiltonian in [2l] 
in a real experiment [5 . Reference [5] discusses such an 
experimental feasibility of this method using cold atoms 
and laser fields. 

Figure [2] illustrates this method for a system composed 
of either one particle or two particles in a one-dimensional 
box performing short range, hard core collisions, |26| with 
the truncated log-oscillator in 27 It reports the proba- 
bility p{Eg) to find the particle at energy Eg during a 
long simulation run. A symplectic integrator was used 
to produce the trajectory of the total system and the 
initial condition was sampled randomly from the shell 
E% t = 5kgT. The numerically computed probability 
(relative frequency) p{E$) is compared to the expected 
Gibbs distribution calculated from [25] according to the 
standard rules of probability theory as 



P(E S 



e- E ^ T n s (E s ) ^ e- E */ T n s (E s ) 

Z(T) f™e- E s/ T Sls(Es)dEs' 

(28) 

where fls(Es) is the density of states of the system. In 
calculating it we neglect the contribution coming from 
the short range interaction, thus obtaining Qs(Es) oc 
E 1 ^ 2 1 , with ri = 1,2 being the number of particles in 

the system. For n = 1 this yields Qs(Es) cx E s 
while for n — 2 we find that Cls(Es) is a constant. The 
agreement between theory and simulations is excellent. 
Further details and discussion can be found in Refs. [5] 
andlTUl 
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FIG. 2: (Color online) Illustration of Method I. Normalized 
probability density function of energy for a system of n par- 
ticles in a ID box performing short ranged collisions, |26| 
with a truncated log-oscillator, |27| of strength fcsT = 15e. 
The total simulation energy is .Etot = 5ksT, the box length 
is L = We Etot/kBT a ~ 1484cr and the log-oscillator cutoff 
length was set to b = a. Black triangles: numerical simula- 
tion with n = 1. Black dots: numerical simulation with n — 2. 
Blue line: Gibbs distribution at temperature ksT = 15e for 
n = 1, as it follows from [28] Red line: Corresponding Gibbs 
distribution at temperature ksT = 15e for n = 2 as it follows 
from |28| This Figure has been provided by Fei Zhan and is 
adapted here from our Ref. 1101 



V. METHOD II 

An alternative method to produce thermostated dy- 
namics is to couple the system to a free particle via a 
logarithmic interaction potential. More explicitly, the 
statement is that the extended Hamiltonian 

(x, p, X, P) = ff s (x, p) + P 2 /2M + k B Tln (| 5 (x, p) - X\/b) 

(29) 

produces thermostated system dynamics, provided the 
(otherwise arbitrary) function g(x, p) induces ergodic dy- 
namics of the total system. 

To demonstrate this, consider the probability 
p(x.,p,X,P) to find the total system at (x, p, X, P). 
Thanks to the ergodic assumption, this is given by the 
microcanonical distribution 

p(x, p, A, P) = S[E tot - #(x, p, A, P)]/n(E tot ) , (30) 

hence: 



P(x,p) 


_ / dXdPS[E tot 


-H s - P 2 /2M - 


-T\n{\g-X\/b)] 




n(£tat) 


(31) 


Making 


the change of variable A' = A 


- ff(x,p), one 


obtains 


irrespective of g\ 


x, p) 




i>(x,p) 


_ J dX'dP5[E tot 


-H s - P 2 /2M 


-k B T\n{\X'\/b)] 




n(E tot ) 


(32) 
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Note that the numerator is the log-oscillator density of 
states fiiog taken at E to t — Hs- Therefore, just as with 
Method I: 



p(x,p) = 



Oiogt^tot - #s(x, p)] e-^W k » T 



Wot) 



Z(T) 



(33) 



The constant temperature equations of motion of this 
second method read: 

x = p/m + k B T[g(x, p) - X^d^gix, p) 
p = -a x C/(x, p) - k B T[g(x, p) - ^"^(x, p) 
X = P/M 

P = k B T[g( X ,p)-X}- 1 

note that for T = the system undergoes constant energy 
dynamics. 

It is important to repeat that thermostated system dy- 
namics are only reached if the global dynamics are er- 
godic. As illustrated below, this requirement is however 
not too restrictive, because we have the freedom to chose 
the function g(x, p). 



Illustration 

To illustrate the method we considered a quartic oscil- 
lator: 



bin, so as to construct a histogram, which, after normal- 
ization gives an approximation to the actual p(Es). A 
similar procedure was followed for the calculation of p(v). 

To begin with we chose g(x,p) = x. Notwithstanding 
the long integration time, the method fails to converge 
to the desired target distributions, see Fig. [3j panel a). 
This means that with the choice of g(x,p) = x, the overall 
dynamics is not sufficiently ergodic to make the system 
sample the canonical ensemble. 

The ergodicity of the overall dynamics can be im- 
proved by choosing a different form for the function 
g(x,p). Panel b) of Fig. [3] reports the result of a dy- 
namical simulation of the same system as in panel a), 
with the same time-step At and simulation time, but 
with g(x,p) = fcx 4 /4, namely we chose g(x,p) as the 
system potential energy. While we found a very good 
agreement between the computed energy probability dis- 
tribution function and the Gibbs distribution, the agree- 
ment between the computed absolute velocity distribu- 
tion and the target Maxwell distribution is still not very 
good. With g(x,p) — sin(fca; 4 /4), see in panel c) of Fig. 
[3j reasonably good agreement between simulation and 
Maxwell distributions was achieved, while the agreement 
between the energy distribution and the Gibbs distribu- 
tion is excellent. Excellent agreement is achieved with 
longer simulation times, see in panel d) of Fig. 3. 



H s = p 2 /2m + kx 4 /4 



(34) 



We simulated the compound system dynamics us- 
ing a symplectic integrator with a time step At = 
10~ 2 b^M/k B T for a total simulation time T = 1.287 x 
10 9 Ai time steps. In our simulations we set fcgT, b and 
M as units of energy, length and mass, respectively. We 
took (xo, Xo,po, Po) = (2,-1,1,-1) as the initial con- 
dition, k = k B Tb~ 4 , and m = M. We computed the 
probability distribution function p(E$) to find the sys- 
tem at energy Eg, and compared it with the target Gibbs 
distribution, [28l The latter reads 



P(E S ) 



e -Es/k B T E -m 



f °° dE s e- E s/kBT Ei 



-1/4 



(35) 



— 1/4 

where the factor E s stems from the density of states 
of the quartic oscillator: J dxdpS[Es — p 2 / 2m + kx 4 / 4} cx 
Eg 1 ^ 4 . We further computed the probability distribution 
function to find the system with a velocity of modulus 
v, and compared it to the target Maxwell distribution, 
reading: 



p{v) 



e -mu 2 /2fc B T 

/o°° e- mv2 / 2kBT 



(36) 



Following Ref. 02 the numerical evaluation of p{Es) 
proceeded by recording the value of Es, once every 100 
time steps. We divided the energy interval [0, E tot ] in 50 
bins, and counted how many times Es was within each 



VI. REMARKS 

As emphasized above, ergodicity of the global dynam- 
ics constitutes the crucial prerequisite for the presented 
methods to work properly. Ergodicity suffices and no 
stronger condition, e.g., the system being mixing, |18j is 
necessary because all that is needed for the system to 
sample the Gibbs distribution is that the compound sys- 
tem samples the microcanonical distribution. It should 
also be mentioned that ergodicity is a sufficient but not 
necessary condition for the methods to work, namely in 
some cases the methods might work even if ergodicity 
does not hold. 

In Method I, whether ergodicity holds depends on the 
specific choice of the interaction energy /i(x, X), which 
must be chosen in any case weak. In Ref. [5] h(x., X) 
was chosen as a hard-core, short range repulsive inter- 
particle potential, [26j and that was sufficient for achiev- 
ing thermostating. In Method II, the ergodicity prop- 
erty depends on the choice of g(x, p), which in turn fixes 
the interaction term k B Tln |g(x, p) — X\. It must be 
emphasized however that our analysis does neither show 
formally nor numerically that the total dynamics are in- 
deed ergodic in the examples presented, but only that, 
loosely speaking, the system appears "ergodic enough" 
for the methods to work. 

Note that in method II the interaction term 
k B Tln \g(x, p) — X\ gives rise to long-range forces. So 
at variance with the implementation of Method I in Ref. 
[5] where the system and the "bath" interacted sporadi- 
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cally through almost instantaneous collisions, in Method 
II they constantly influence each other, due to the long 
range force. 

We have shown how different choices of <?(x, p) can 
result in different ergodic properties of Method II. An 
important subject for further studies would be to derive 
a set of criteria for appropriately choosing g(x, p), given 
the properties of the system, as encoded in its Hamilto- 
nian ffg(x, p). 

Besides choosing /i(x, X) or g(x, p), the ergodicity 
of both methods can be improved also by substituting 
the log-oscillator with a multi-dimensional log-oscillator, 
which will add more degrees of freedom to the whole sys- 
tem, see in Appendix. 

In implementing Method II, we have replaced the loga- 
rithmic potential with the same truncated potential, 27 
used for Method I. Therefore, just as with Method I, this 
truncation can lead to deviations to the target Maxwell 
distribution when the number of particles in the system 
increases. An interesting line for future studies would 
then be to put forward implementations that avoid the 
truncation and treat the singularity by some other means, 
which might allow for applying the methods to large sys- 
tems as well. 



VII. CONCLUSIONS 
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FIG. 3: Illustration of Method II. Each panel reports the 
analytically and numerically computed probability distribu- 
tion for system energy Es (rescaled by the total fixed energy 
Etot) and system speed v (rescaled by the maximal speed 
^/2E to t/m), for various choices of g(x,p). Panel d) 
has the same g(x,p) as panel c), but for a longer simulation 
time. 



With this study we presented two Hamiltonian schemes 
which allow a system H$ to sample a canonical Gibbs dis- 
tribution. This being so, the method of thermostating 
is achieved here in a deterministic time-reversal invari- 
ant and symplectic manner. Both schemes rest upon the 
spectacular thermodynamic property of logarithmic os- 
cillators of having an infinite heat capacity. Hence, in our 
methods a single log-oscillator substitutes an infinite heat 
bath coupled weakly to the system. With our Method I 
we couple the system weakly to a log-oscillator where the 
absolute temperature T denotes the strength of the loga- 
rithmic potential. In Method II we consider a composite 
system of H$ and a free particle which is coupled with 
a long range log-interaction of strength T to the system 
of interest H$- Note that Gibbs thermalization occurs 
here independent of the interaction-strength T, being ei- 
ther strong (large T) or weak (small T). A prominent 
property inherent to both schemes is that these are man- 
ifestly Hamiltonian QJ. Also, at variance with the Nose 
Hamiltonian, [l] our Hamilton functions possess standard 
(i.e. coordinate-independent) kinetic energy contribu- 
tions. This fact in turn allows not only an implementa- 
tion with numerical means but as well a physical realiza- 
tion. This advantage should be contrasted nevertheless 
with the limitation that both methods inherit from per- 
forming a truncation of the logarithmic potential as in 
|27[ which, as thoroughly emphasized in our previous ac- 
counts [3 [6] , limits an efficient thermostating to systems 
with a small number of degrees of freedom. Notably, the 
investigation of such nano-scale systems is in the lime- 
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light of present day research activities 
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Appendix A: Appendix. / dimensional 
log-oscillators 

Consider a / dimensional log-oscillator: 

tf log (X,P) = |^+/fc B Tln^, (Al) 

Where X = (X u X 2 , -X f ), P = (P u P 2 , ...P f ). For the 
phase volume &\ og (E) = j H<E dXdP one obtains: 



(8TT 2 b 2 Mk B T/f) 

r(/ + i) 



//2 



E/k B T 



(A2) 



where Y denotes the Gamma function. Therefore, the 
density of states is exponential in E/k B T; reading 



d<t> log (E) _ (8ir 2 b 2 Mk B T/f) 



//2 



dE 



k B TY(f + l) 



?E/k B T 



(A3) 



Consequently, the methods presented above can also be 
implemented with an / dimensional oscillators replacing 
the 1 dimensional oscillator: In this case Method I be- 
comes: 



x = p/m, 

p = -9 x C/(x)-9 x /i(x,X) 
X = P/M 

P = -[k B T/X 2 ]X~ dxh(x,X) 



and Method II becomes 

x = p/m + [k B T/( E - X) 2 ] J2 k (g k - X k )8 p g k 
P = -dJJ - [k B T/(g - X) 2 ] J2 k (g k - X k )d x g k 
X = P/M 

P=[fc s T/(g-X) 2 ](g-X) 

where g, a short notation for g(x, p) = 
(<?i(x, p), . . . , <?/(x, p)), is an /-dimensional field. 
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